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of the second Aeroelastic Prediction Workshop for the NASA Benchmark Supercritical Wing configuration. 
The computational results are obtained using FUN3D, an unstructured grid Reynolds-Averaged Navier-Stokes 
solver developed at the NASA Langley Research Center. The analysis results focus on understanding the dip 
in the transonic flutter boundary at a single Mach number (0.74), exploring an angle of attack range of —1° to 
8° and dynamic pressures from wind off to beyond flutter onset. The rigid analysis results are examined for 
insights into the behavior of the aeroelastic system. Both static and dynamic aeroelastic simulation results are 
also examined. 
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I. Introduction 


eometric simplicity first drew us to the Benchmark Supercritical Wing (BSCW) as an analysis configuration for 
Gi. Aeroelastic Prediction Workshop (AePW). With an outer mold line machined from a block of aluminum, 
the wing itself is one of the most rigid aeroelastic models ever tested in the Transonic Dynamics Tunnel (TDT); the 
flexibility was provided principally by the TDT Pitch and Plunge Apparatus (PAPA). The PAPA is a two degree-of- 
freedom (pitch and plunge) sidewall mount system constructed of bars arranged to produce the desired pitch and plunge 
mode stiffness values, and masses positioned to decouple the wind off modal mass properties. This mount system 
was engineered by Farmer,! who reported experimental results for a supercritical wing showing similar transonic 
aeroelastic characteristics to the computational results presented in the current paper. 

The simplicity of the PAPA structure has allowed us to study the flow field evolution as Mach, angle of attack and 
excitation level are varied, and the resulting effects on the aeroelastic deformation (static aeroelastic deflection) and 
stability (flutter). This paper and its companion publication? endeavor to present a sufficient summary of results to 
show the influences of shock development and separation onset on the aeroelastic characteristics of the BSCW wing. 
The goal of the work is to develop and understand the bucket in the transonic flutter boundary for this configuration. 
Historically, the location (in terms of Mach number) and depth (in terms of dynamic pressure or velocity) of the 
transonic dip have been shown to be functions of the angle of attack. There are many configurations where this 
behavior has been observed in the technical literature, including both computational and experimental examples. A 
specific example for the NACA 64 A-010 airfoil is reproduced from reference 3 in Figure 1. 
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Figure 1. Effect of angle of attack on NACA 64 A-010 flutter boundaries; Source: Reference 3, Figure 7. 


While previous publications associated with AePW-2 concentrated on the overall summary of the computational 
results generated for the workshop,’ in this paper, we will be reporting on the predictions from a single code, FUN3D. 
We focus principally on the influences of the separated flow onset, induced in this study through three effects: angle 
of attack, dynamic pressure and initial excitation. Only one Mach number is examined—Mach 0.74. Angles of attack 
from -1° to 8° are examined, with the primary investigations at angles of attack of 0°, 3.15°, 5° and 6.5°. The dynamic 
pressure ranges examined for each angle of attack begin at the rigid condition and range to above the predicted flutter 
onset condition. Distributions of steady pressure coefficient for the static rigid BSCW configuration are compared 
with experimental information over the angle of attack range where that information is available, which ends at 5°. 
The flutter onset prediction at 0° is compared with the corresponding result from experiment. Unfortunately, the 
experimental flutter boundary at nonzero angles of attack is not available for a similar comparison at this Mach number. 
In fact, very little quantitative information is available for any of the conditions where the simulations indicate that 
separated flow effects dominate the stability. 

The paper is organized into five major sections. Background material on the configuration and its uses in the 
Aeroelastic Prediction Workshops are presented first. The analysis processes are then discussed, followed by results 
from each of the three analysis processes used here: rigid steady, static aeroelastic and dynamic aeroelastic analyses. 
The dynamic aeroelastic results include effects of matched point conditions, single degree of freedom representa- 
tion, perturbation level and spatial resolution. The final two sections contain discussions of continuing analyses and 
conclusions. 
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II. Background 


A. The BSCW Model and Test Cases 


The BSCW model, shown in Figure 2, has a simple, rectangular, 16- x 32-inch wing planform, with a NASA SC(2)- 
0414 airfoil. It was first tested in the TDT in 1991.° For this test, the wing was mounted on the TDT Pitch and Plunge 
Apparatus (PAPA) to obtain the flutter boundary at various Mach numbers and angles of attack for this two-degree of 
freedom system. In 2000, the wing was tested again, this time on an Oscillating Turntable (OTT).° The purpose of the 
OTT test was to measure aerodynamic response during sinusoidal (forced) pitch oscillation of the wing. 

For both the OTT and PAPA tests, the model was mounted to a large splitter plate, sufficiently offset from the wind- 
tunnel wall (40 inches) to (1) place the wing closer to the tunnel centerline and (2) be outside the tunnel wall boundary 
layer.’ The wing was designed to be rigid, with the following structural frequencies for the combined installed wing 
and OTT mounting system: 24.1 Hz (spanwise first bending mode), 27.0 Hz (in-plane first bending mode), and 79.9 
Hz (first torsion mode). When installed on the PAPA mount, the combined system frequencies were 3.33 Hz for 
the plunge mode and 5.20 Hz for the pitch mode.* The plunge and pitch modes are the only modes included in the 
aeroelastic analyses of the PAPA-mounted configuration to be presented in this paper. 

For instrumentation, the model has pressure ports in chordwise rows at the 60% and 95% span locations. For 
the BSCW/PAPA test, both rows were fully populated with unsteady in situ pressure transducers. The quantitative 
information obtained consists of unsteady data at flutter points and averaged data on a rigidized apparatus at the 
flutter conditions. For the BSCW/OTT test, only the inboard row at 60% span was populated with transducers. The 
quantitative information for the OTT test consists of unsteady pressure data and accelerometer data for forced pitch 
oscillations and for the unforced system at stabilized flow conditions. 


¥ 
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(a) Photograph of the BSCW model mounted on the OTT in 
the TDT. 


(c) Cross-sectional view of the SC(2)-0414 airfoil, with BSCW in- 
strumentation. 


Figure 2. BSCW model. 


B. The First Aeroelastic Prediction Workshop 


The first AIAA Aeroelastic Prediction Workshop (AePW-1) took place in conjunction with the AIAA Structural Dy- 
namics and Materials (SDM) conference in April 2012.° This workshop utilized three configurations, including the 
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BSCW. There were three BSCW test cases for this workshop: a rigid steady analysis and two forced oscillation anal- 
yses, all at Mach 0.85, 5° angle of attack. The submitted simulation results for the BSCW contained the most widely 
varying results of all of the AePW-1 configurations.!° These variations were observed in the steady results as well as 
in the frequency response functions. Predictions of the upper surface shock location spanned more than 20 percent of 
the chord.'! Because of these wide variations, further studies were initiated to understand the underlying causes. 


C. The Second Aeroelastic Prediction Workshop 


The second AIAA Aeroelastic Prediction Workshop (AePW-2) took place in conjunction with the AIAA SciTech 2016 
Conference.* The computational aeroelasticity community was challenged to analyze the BSCW configuration at three 
conditions (M=0.7 at 3° angle of attack on the OTT; M=0.74 at 0° angle of attack on the PAPA; and M=0.85 at 5° angle 
of attack on the PAPA) and to present their results at the workshop. The experimental data from the OTT test indicated 
that the BSCW exhibited a strong shock and boundary-layer-induced separated flow at a moderate angle of attack at 
transonic conditions.!* The computations of the transonic flow, in conjunction with the flutter boundary predictions, 
were therefore the focus of the workshop.'? The summary of flutter onset predictions from AePW-2 is shown in 
Figure 3. The predictions at Mach 0.74, 0° angle of attack have considerably less scatter than the prediction at Mach 
0.85, 5° angle of attack. The lower Mach number, lower angle of attack predictions also surround the experimental 
flutter onset condition, 169 psf. 
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Figure 3. Flutter boundary results from AePW-2. 


D. Grids 


For the present study, two of the unstructured grids that are detailed in reference 2 were utilized: the designated 
“medium” grid (9 million nodes) and the designated “Grid D” (35 million nodes). In the current paper, Grid D is 
referred to as the extra fine grid. The grid distributions for both the wing surface and the plane of symmetry are 
presented in Figure 4. 
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(a) Medium grid; Grid B; 9 million nodes. (b) Extra fine grid, Grid D; 35 million nodes. 


Figure 4. Medium and extra fine grids. 


II. Analysis Processes 


Solutions to the Reynolds-Averaged Navier-Stokes (RANS) equations were computed using the NASA Langley- 
developed computational fluid dynamics (CFD) software FUN3D.!* The flow solutions were produced with turbulence 
closure obtained using the Spalart-Allmaras (SA) one-equation model.!>-!° Inviscid fluxes were computed using the 
Roe scheme,!” primarily without a flux limiter. The solutions were produced by running three distinct analyses, 
differentiated by the treatment of the structure, the initial condition and the temporal parameters. The rigid, static 
aeroelastic and dynamic aeroelastic solution processes are described below. 


A. Rigid Analysis 


The initial solution for each angle of attack at Mach 0.74 was generated for a rigid wing. In general, these solutions 
were run ina “steady” mode, in which the solutions computed were not time-accurate and the structure was represented 
as a fixed position outer mold line. The position (angle) was different for each angle of attack and analyzed in an 
individual simulation. Time integration was accomplished by an Euler implicit backwards difference scheme, with 
local time stepping to accelerate convergence. Although some of the cases were not asymptotically steady, the rigid 
steady simulation was still the first step performed. The cases that were not asymptotically steady occurred at the 
higher angles of attack and will be identified as a subset of the results presented. These cases did not converge to a 
fixed pressure distribution and should be run in a time-accurate manner where the range of the physical solution can 
be established. Results of time-accurate rigid simulations are not presented in the current paper. 


B. Unsteady Aeroelastic Analyses 


For aeroelastic analyses, the flow equations were solved in a time-accurate mode!* coupled with a structural solver, 
allowing the grid to deform. For this work, this coupling was performed internal to the FUN3D code!? using a 
modal representation of the structural dynamics. This modal solver was formulated and implemented in FUN3D in a 
manner similar to other NASA Langley aeroelastic codes (CAP-TSD”” and CFL3D*!). For the BSCW computations 
presented here, the structural modes were obtained via normal modes analysis (solution 103) with the Finite Element 
Model (FEM) solver in MSC Nastran?™ 2? The modes were then interpolated to the surface mesh using the method 
developed by Samareh.*? The BSCW FEM is described in reference 13. 

To solve the unsteady flow equations, the FUN3D solver employs a dual-time-stepping method, which is widely 
used in CFD. This method involves adding a pseudo-time derivative of the conserved variables to the physical time 
derivative that appears in the Navier-Stokes equations. Subiterations are run within each physical (global) time step 
to converge the flow solution. Convergence here refers to the process where the pseudo-time derivative vanishes as 
the subiterations proceed within each physical time step. At the end of the iterative process, a solution to the original 
unsteady Navier-Stokes equations is obtained for that physical time step. For an infinitely large physical time step, the 
dual-time solver becomes identical to the steady-state solver, though of course, all time accuracy is lost. 

Aeroelastic analysis requires a grid deformation capability. The grid deformation in FUN3D is treated as a linear 
elasticity problem. In this approach, the grid points near the body can move significantly, while the points farther away 


5 of 24 


American Institute of Aeronautics and Astronautics 


may not move at all. For the aeroelastic simulations in this paper, a predictor-corrector coupling method between the 
flow solution and the structural solution has been used. Within each physical time step, the flow field solution is 
converged using the subiterations and subsequently coupled with the structural equations. This coupling occurs only 
once for each physical time step. The coupling is described further and explored in detail in reference 2. 


I. Static Aeroelastic Analysis 


The static aeroelastic simulation procedure begins from the last point in the rigid steady solution at identical Mach 
number and unloaded (rigid) angle of attack. Static aeroelastic simulations are obtained by restarting the CFD analysis 
from this rigid steady solution. A high value of structural damping (0.9999) is used to allow the structure to find 
the equilibrium position with respect to the mean flow. The fluid-structure coupling provides a step input to the 
system, initiating the migration from the rigid solution to the coupled static aeroelastic solution. In general, the static 
aeroelastic analysis is performed using larger time steps and fewer subiterations than have been shown to be required 
for a time-accurate temporally-converged solution. This issue is briefly explored later in this paper. 


2. Dynamic Aeroelastic Analysis 


Dynamic (flutter) analysis is performed in a multistep process, with the dynamic solution started from the static 
aeroelastic solution at corresponding Mach number, angle of attack and dynamic pressure, and setting the structural 
damping value to zero. An excitation is provided to the system by adding an initial condition to the generalized 
(modal) velocities. For all results presented in the current paper, identical values of initial generalized velocity are 
applied simultaneously to both the plunge and pitch modes. The effects of the initial excitation value on the flutter 
solution are discussed later in this paper. 

The wing response, in the form of the time varying pitch angle, is computed and a logarithmic decrement method 
is used to calculate the damping ratio. For a stable solution, the damping ratio is greater than zero, and for an unstable 
solution, the damping ratio is less than zero. To determine the flutter onset dynamic pressure, linear interpolation is 
performed to find the zero damping case using the two cases on either side of the stability boundary (that is, using 
the two dynamic pressures where the sign of the damping changes in between them). The flutter frequency is then 
calculated using the same interpolation coefficients. 

Each simulation performed has a different number of time steps and each case has different transient behavior. 
These two factors mean that the uncertainty associated with each calculated damping value will be different. In 
the current paper, the uncertainty is captured by standard error values, which are presented along with the damping 
values in some cases. The standard error of the estimate of the damping value is computed based on the logarithmic 
decrement process information. The logarithmic decrement technique, as applied, performs a linear fit to the logarithm 
of the extreme values of the oscillations (that is, the peaks and valleys of the cycles). The damping is computed as in 
equation 1, where m is the slope of the linear estimate and q@ is the estimate of the frequency in radians. 


m 
c=” (1) 


The standard error estimate of the slope from a linear regression model is calculated as in equation 2, where x; 
are the values of the independent variable (time), and y; are the corresponding values of the dependent variable (the 
logarithm of the extreme values of the pitch angle). ¥; are the estimated values of the dependent variable computed 
using the linear regression model (that is, using the slope and intercept computed from the original data pairs (x;,y;)). 
N is the number of degrees of freedom (the number of peaks and valleys used in the logarithmic decrement process). 
Note that the number of degrees of freedom is reduced by 2, since 2 parameters (slope and intercept) were estimated 
to formulate the sample estimates, yj. 
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The standard error estimate of the damping value is then computed from the standard error estimate of the slope, 
assigning no error in the calculation of the frequency, equation 3. 


SE¢ = SEn/o 3) 
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IV. Simulation Results 


Rigid, static aeroelastic and dynamic aeroelastic analyses were conducted at Mach 0.74 to investigate the influ- 
ence of separation onset. The angle of attack was was varied from —1° to 8° and different excitation levels were 
applied to the system. Multiple dynamic pressures, each analyzed through an independent simulation, were analyzed 
at combinations of angle of attack and excitation level. 


A. Rigid Steady Simulations 


Results from rigid simulations are presented in this section. These results will show that the flow field characteristics 
change when the airfoil exceeds 6°. This change is shown through surface and volumetric plots (Figure 5), integrated 
lift coefficients (Figure 6), and pressure field slices (Figures 7 and 8). 

In the rigid steady simulations, as the angle of attack increases, the region of flow separation increases and the 
shock structure moves forward as shown in Figure 5. The rigid steady simulations produce converged results for the 
simulations below the critical angle of attack of 6°; above this angle, the solution meanders rather than converges. This 
meandering behavior is illustrated in Figure 6, which shows the development of the lift coefficient for different angles 
of attack as the iterations proceed. In Figure 6, the iteration histories corresponding to the highest two angles of attack 
(Gp = 6.5° and 8°) do not converge to fixed values. This nonconvergent behavior corresponds to conditions where 
the flow field is extensively separated and which are outside the range where we think that URANS analysis produces 
credible results. As previously noted, experimental information does not exist for these higher angles of attack, so this 
latter assessment is admittedly made without being able to benchmark the computational results. 


(a) 0° 


(d) 6.5° (e) 8.0° 


Figure 5. Skin friction surface plots from rigid analysis and Mach contour plots at splitter plate; Snapshots at end of 
simulation. 


Figure 7 shows the pressure coefficient versus nondimensional chord location at the 60% span station on both the 
upper and lower surfaces for these same angles of attack. In these plots, and all pressure coefficient plots presented, the 
vertical axis is inverted. As shown in the figure, the shock moves aft with increasing angle of attack for angles below 
6.5° and then reverses direction as the angle of attack is further increased. The upper surface shock forms between 
1.15° and 2.4°; the shock continues to gain strength and move aft until 6°. For higher angles, the shock then moves 
forward with increasing angle of attack, still strengthening. 

Aft of the shock, a small region of separation is indicated for angles of attack of 5° and higher, as illustrated in the 
skin friction plots in Figure 5. In these plots, skin friction is shown on the surface of the wing, with orange indicating 
positive skin friction (attached flow) and blue indicating negative skin friction (separated flow). On the upper surface, 
the plots for 3.15° and below indicate no separated flow region. For higher angles, the blue regions indicate the extent 
of the separated flow region. In these pictures, a slice through the flow field is also shown at the intersection between 
the splitter plate and the wing root. The contours on this slice show the local Mach number, emphasizing the upper 
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Figure 6. Rigid simulation as iterations proceed; Medium grid. 


surface shock location through an abrupt color change, generally from red to purple. The shock extends across the 
wing, with a small region of separation behind the shock for cases of 5° and above. As the angle of attack increases, the 
figures show that near the splitter plate, the shock progresses toward the leading edge and is strengthened, indicated 
by the progression of the color shown in the shock from red to orange to yellow. Near the splitter plate, there is 
no observed reversal of the direction of shock travel, but the shock appears to weaken near the splitter plate for cases 
above 6.5°. This is in contrast to the pressure data shown in Figure 7 at the 60% span station. In addition, the separated 
flow region shown near the splitter plate in Figure 5 is very small in extent in all directions. Starting at 6°, however, 
there is a larger region of separated flow, shown by the growth of the blue region on the wing as the angle of attack 
increases. 

The plots in Figure 5 also indicate the effect of the tip vortex. The grey streamers in the flow are streaklines, as 
are the black lines on the slice at the splitter plate. Study of the streaklines shows that the flow over the aft portion of 
the wing upper surface flows inboard, as well as streamwise. On the lower surface, the streaklines show that the flow 
moves outboard as well as streamwise. Although the influence of the tip vortex has not been examined in detail in the 
current work, limited studies were conducted at @ = 5° on a two-dimensional airfoil section. This work is incomplete 
at this point in time and the preliminary results are not presented here, but it was evident that the two-dimensional 
results were qualitatively different from the three-dimensional results. 
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(a) Upper surface (b) Lower surface 
Figure 7. Final time point snapshots of pressure coefficients at 60% span station; Medium grid; Rigid steady analysis. 
Comparison of the computational results with experimental results is shown in Figure 8 for selected angles of 
attack where experimental quantitative information is available. The computational results were taken as snapshots 


at the last time point from each simulation and are shown by the lines; experimental information is represented by 
the symbols. The grey shaded regions show the range of the experimental pressure data at the flutter point. The 
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computations predict the upper surface shock to be stronger and further aft than indicated by the experimental data. 
This trend is consistent with the trends observed in the AePW data sets.+!' Although the previous work has identified 
the lower surface aft region as being mispredicted, this trend isn’t so much in evidence in the current calculations. 
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Figure 8. Pressure coefficients at 60% span station; Rigid steady case. Symbols: Experiment. Colored Lines: Analysis. The 
grey region shows the range for the flutter information from the experiment; Symbols show the mean. 


The lift and pitching moment coefficients are shown as functions of angle of attack in Figure 9 for the medium 
grid. Both the lift and the pitching moment coefficient show linear dependence on the angle of attack until 6°, where 
dramatic changes in both curves are seen and stall is indicated. 

The BSCW airfoil is a supercritical airfoil (SC(2)-0414), so it was designed with the goal of being shock-free at 
a given lift coefficient; the design lift coefficient for this airfoil is 0.4. Note that at the Mach number of interest for 
this study, 0.74, a lift coefficient of 0.4 is produced by an angle of attack near 2°. The pressure coefficient plots in 
Figure 7 show that the shock forms between 1.15° and 2.4°, correlating well with the designed behavior of the airfoil. 
The coefficient values presented are the final time point of each simulation, not the mean value nor the range of the 
coefficients. Statistical values should be obtained from time-accurate simulations to better capture the behavior that 
the method is predicting. 

The aerodynamic center and the center of pressure were computed from the coefficient data. The center of pressure 
is shown in Figure 9c, measured positive aft with the origin at the wing leading edge. The wing chord, 16 inches, was 
used as the reference length. Thus, the trailing edge location is at a value of | in the figure. The center of pressure 
was calculated from the rigid data, using equation 4 after applying the parallel axis theorem to reference the pitching 
moment coefficient to the wing leading edge. (The pitching moment coefficient is defined as positive leading edge 
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Figure 9. Steady rigid simulation results; Medium grid. 


up. The lift coefficient is defined as positive upward. The center of pressure location as defined in equation 4 is thus 
positive aft from the leading edge.) 


= Cin 
Xcp = Cc (4) 
The center of pressure is aft of the trailing edge of the airfoil for negative angles of attack, as shown in Figure 9c. At 
0° angle of attack, the center of pressure is aft of the midchord and moves forward as the angle of attack increases. The 
forward-most position of the center of pressure occurs at 6° and is at 34.6% chord behind the leading edge. Further 
increases in angle of attack result in moderate migration of the center of pressure aftward. The center of pressure 
migration is similar to the movement measured and reported in the earliest of airfoil tables for asymmetric airfoils, 
particularly paralleling those with lower surface cusp regions in the aft portion of the airfoil.7* 
The aerodynamic center was also calculated from the rigid steady data sets using equation 5. 


C ma 
X iC SS == — 5 
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The results are shown in Figure 9d. The pitching moment coefficients were calculated relative to the leading edge, 
thus the location values given are also with respect to the leading edge. In the range where Cz and C,, are linear 
functions of the angle of attack, the aerodynamic center location is constant (21.5% chord aft of the leading edge). 
In the region where the behaviors of the coefficients are nonlinear, the aerodynamic center location shifts erratically. 
Calculations at a refined set of angles of attack in the nonlinear response range would be required to better characterize 
this movement. The aft shift in the aerodynamic center is produced by separating flow, which will be shown in a later 
section of this paper to have a destabilizing effect on the system. 
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B. Static Aeroelastic Results 


Static aeroelastic results were generated using the CFD analysis process outlined previously. The static aeroelastic 
solutions did not develop monotonically and in some cases looked converged for a large number of iterations before 
migrating to a final equilibrium position. While most of these static aeroelastic solutions were performed with a time 
step and subiterations that were insufficient to characterize each time step as time-accurate (i.e., the time step was too 
large and the subiterations too few to declare that the time history points had physical significance), a refined time step 
solution for & = 5°, Gg = 152 psf will be examined later in this section. 

Static aeroelastic solutions achieve different equilibrium angles of attack based on the rigid (wind-off) angle of 
attack, the dynamic pressure, the Mach number and the Reynolds number. The latter two parameters are fixed values 
for this study. Increasing dynamic pressure induces additional twist as it amplifies the influence of the pressure relative 
to the structural stiffness. This effect is captured in Figure 10, which shows the medium grid predicted aeroelastic twist 
increment as a function of dynamic pressure for rigid angles of attack from 0° to 5°. The figure shows results from the 
CFD simulations using circular symbols; square symbols show results from simplified analyses that will be discussed 
later in this section. In static aeroelastic solutions for rigid angles of attack less than or equal to 1.15°, dynamic 
pressure induces negative (nose down) aeroelastic twist, reducing the angle of attack. For angles of attack greater than 
1.15°, the induced aeroelastic twist angle is positive (nose up) and is linear with dynamic pressure until the total angle 
of attack exceeds 6.5°. Beyond this condition, additional increase in either dynamic pressure or angle of attack results 
in increasing the separation region and the shock structure moving forward, producing a corresponding decrease in the 
total angle of attack. The data sets presented in Figure 10 for a = 4° (green line) and a@& = 5° (grey line) exemplify 
the influence of separated flow onset, as nonlinear behavior is observed beginning near 200 psf for %& =4° and 120 psf 
for Q% = 5°. This will be examined subsequently for the o&% =5° case. Note that the change in induced twist direction 
occurs at the same angle of attack where the upper surface shock forms. That is, for shock-free cases, the induced twist 
angle is leading edge down; for cases after the upper surface shock has developed, the induced twist angle is leading 
edge up. 

The influence of the separating flow is best observed for the G@ = 5° case, which is shown in Figure 10 and 
highlighted in Figure 11. In time history plots of the developing static aeroelastic solutions, there are temporary 
plateaus in the aeroelastic twist angles before the solution transitions to another equilibrium angle. These temporary 
plateau values are shown in Figure 11 by the black-filled circles. The grey-filled circles show the final values obtained 
via the medium grid. The difference in the early versus final solutions is pronounced for dynamic pressures greater 
than 100 psf, where the computed twist angle deviates considerably, creating a new slope of the line. The quantitative 
information shown in this figure also includes simulations using the extra fine grid. The extra fine grid results match 
the early equilibrium (plateau) values from the medium grid solution. It is currently unclear whether the extra fine 
grid results require additional simulation time and will eventually achieve the reduced angle equilibrium condition or 
whether the medium grid results produce misleading trends at the higher dynamic pressures. Although extra fine grid 
results were generated with and without a flux limiter, no significant differences were observed in the static aeroelastic 
solutions. 
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Figure 10. Static aeroelastic twist angles; Medium grid; CFD compared to quasi-steady equation. 
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Figure 11. Static aeroelastic twist angles, 5°. 


To better understand the transition behavior observed, several simulations with refined temporal parameters were 
performed. Figure 12 shows results for analyses performed at OQ} = 5°, 152 psf. The result labeled “Coarse Temporal 
Parameters” corresponds to the nondimensional time step size (DT) = 50, dimensional time step size (dt) = 0.0082 
seconds, with ten subiterations and Courant-Freidrichs-Lewy (CFL) numbers for mean flow equations = 25 and CFL 
numbers for turbulence model equations = 50. These parameters were determined by performing a convergence 
efficiency study at the AePW-2 Case 2 condition (Mach 0.74, % = 0°) and correspond to the parameters standardly 
used for static aeroelastic analyses in the current paper. The result labeled “Fine Temporal Parameters” corresponds 
to a simulation performed with the parameters used in the dynamic simulations (DT = 1.12, dt = 0.0002 seconds, 200 
subiterations maximum, temporal error convergence of 0.10). The two additional simulation results shown are for 
intermediate choice of time step size and reduced CFL numbers. 

In all four cases, the simulations find an initial static twist angle but then migrate to a smaller angle. While the 
smaller twist angle resembles a stabilized position for a considerable number of iterations in most cases, the continued 
simulations all show at least one oscillation back toward the original preferred position. For example, the coarse 
parameter iterations have a hump between 5000 and 5500 iterations where the twist angle changes from 1.4° to 1.5° 
and then back again. This suggests that one shock location is slightly more stable than the other. The more stable shock 
location is the forward location, with the wing at a lower angle of attack. The fine time step result shows decaying 
sinusoidal oscillations with a peak-to-peak range of approximately 0.1° and a down-trending mean value. It is unclear 
whether or not the mean value will eventually reach the equilibrium position indicated by the larger time step values. 
Examination of the flow field details makes it clear that the higher angle static aeroelastic simulations correspond to 
forward shock locations and larger separated flow regions, however, it is unclear what aspect of the solution triggers 
the transition from one flow state or shock location to the other. 
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Figure 12. Static aeroelastic twist angles, % = 5°, 152 psf; Medium grid; Varying temporal solution parameters. 
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The dynamic pressure influence on the shock position is illustrated for the case of the medium grid analyses at 
QO = 5° in Figure 13. As the increasing dynamic pressure initially twists the airfoil, the shock moves aft. However, 
once the airfoil twists such that the total angle of attack is greater than approximately 6.5°, the shock moves forward. 
This twist exceedance and forward motion is shown in Figure 13 to correspond to cases where the dynamic pressure 
is greater than 100 psf, which matches the condition where the transitory equilibrium behavior begins. 


(a) 60% span. (b) Upper surface shock region. 


Figure 13. Pressure coefficient distributions from static aeroelastic simulations for various dynamic pressure, % = 5°. 


Single degree of freedom estimates of the static aeroelastic twist were also generated based on the rigid static 
computations. These estimates are useful for trying to understand which cases require analysis using the coupled 
simulation capabilities and which can be derived from the much-simplified steady rigid procedure. In the case of an 
aeroelastic system, the total angle of attack is not the rigid angle of attack, but rather is the rigid (wind-off) angle of 
attack plus the aeroelastically induced twist angle, equation 6. The loads are found by using the linear fit estimates of 
the aerodynamic derivatives. The flexible static response using the linear method is calculated by setting the aerody- 
namic moment equal to the structural restorative moment. This results in an expression for the aeroelastic twist angle 
given in equation 7, where & is the rigid angle of attack plus a bias angle found as the x-intercept of the plot of pitch- 
ing moment coefficient versus angle of attack. The comparison results are shown in Figures 10 and 11. Examining 
the results, it is only post-stall that the coupled CFD static aeroelastic characteristics deviate significantly from those 
derived from the single degree of freedom rigid simulations. 


Aoral = Oy + 0 (6) 
Z 1 
q—S*c*Cma—1 


C. Dynamic Aeroelastic Simulations 


Dynamic aeroelastic simulation results are produced by applying an initial perturbation of generalized velocity to the 
coupled system equations, as described in Section B. The coupled system is observed to respond in an oscillatory 
manner, and the time histories are analyzed. For the BSCW at Mach 0.74, there are different flow field scenarios in the 
dynamic simulations: attached flow, separated flow and separating/attaching flow. Summary plots of the medium grid 
flutter results are shown in Figure 14 for four angles of attack (0°, 3.15°, 5° and 6.5°) and two values of generalized 
velocity excitation (gvelg =0.5 and 5) . The root locus plot, Figure 14a, shows the progression of the pitch eigenvalue 
as dynamic pressure increases. Figures 14b and c show damping and frequency values, respectively, versus dynamic 
pressure for the medium grid simulations. Time histories of the aeroelastic pitch angle increment, 0, for flutter onset 
conditions are shown in Figure 15 for rigid angles of attack of 0°, 3.15° and 5°. Time histories of pitching moment 
coefficient and @ for the same angles and also 6.5° are shown in Figure 16. The pitching moment coefficient time 
histories (blue lines) emphasize the change in the physics that occurs for the Q% = 6.5° case, shown by the nonlinear 
responses. The aeroelastic pitch angle increments vary sinusoidally and smoothly throughout, as shown by the green 
lines in Figure 16. 
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Figure 14. Flutter results; Medium grid 
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Figure 15. Time histories of analyzed conditions above and below the flutter conditions; Medium grid. 


For the 0° angle of attack case, the flow is attached without an upper surface shock. Here, the static aeroelastic 
solution (see Figure 10) induces approximately 1° nose down twist, reducing the total effective angle of attack. All 
of the time histories (see Figures 15a and 16a) appear perfectly sinusoidal. The resulting flutter onset as seen in 
Figure 14 occurs at a dynamic pressure of 152 psf, which is 10% below the experimentally measured onset dynamic 
pressures of 169 psf. The pitch mode frequency is seen to migrate smoothly and monotonically from the rigid pitch 
mode frequency at 5.2 Hz, down to 4.2 Hz at flutter onset. The damping trend is similarly smooth, with the maximum 
damping calculated at 100 psf. The flutter onset results are nearly identical for both values of perturbation velocity, 
gvelg, shown in Figure 14. 

The behavior of the pressure field for different parts of the oscillation cycle is shown in Figure 17. For this case, 
there is no upper surface shock and no separated flow, but a shock does appear on the lower surface at the extreme 
negative angles of attack. Information is presented at three dynamic pressures in Figure 17: below flutter onset (subplot 
a, g = 135 psf); the closest case to flutter onset (subplot b, g = 152 psf) and post-flutter (subplot c, g = 169 psf). In 
each subplot, there are four graphics. The top graphic shows the time history of the aeroelastic twist angle increment, 
6, with six specific time points identified at the oscillation peaks and valleys. The second graphic shows the snapshots 
of the pressure coefficient distribution at 60% span as a function of chord. Snapshots of the upper surface pressure 
distribution are only shown at the time points identified in the top figure with colored symbols. The line colors in 
the second graphic correspond to the symbol colors in the first graphic. Snapshots of the lower surface pressure 
distribution are shown for all time points with the same color-coordination as the upper surface snapshots. Snapshots 
corresponding to time points that aren’t specially identified by symbols in the top graphic are colored brown. 

The dominant feature of the snapshots is the brown region that shows the range of the lower surface pressures. 
Each valley identified in the oscillation, corresponding to the lowest aeroelastic twist angle and thus the lowest angle 
of attack, corresponds to a lower surface shock at the aftmost position on the chord. The shock position repeats nearly 
identically for the three identified valley points. The third graphic in each subfigure shows a zoomed region of the time 
history of the aeroelastic twist angle, focused on the time range with the points previously identified. These identified 
time points are aligned with the time axis in the fourth graphic. The fourth graphic is a contour plot of the pressure 
coefficient on the upper surface, shown as a function of chord (vertical axis) and time (horizontal axis). For the 0° 
case, the upper surface behavior is very benign, and the contour plots show little. 

Attached flow with an upper surface shock is seen in the case of @& = 3.15°. Here, the static aeroelastic solution 
shows an increased angle of attack as dynamic pressure increases in Figure 10. Static aeroelastic solutions were only 
performed at dynamic pressures up to 169 psf. At this maximum condition, the total angle of attack was less than 6.5°. 
The dynamic solutions, using gvelp = 5.0 produced a flutter onset condition nearly identical to that produced at 0°, as 
shown in Figure 14. Time histories of integrated quantities are all nearly perfect sinusoids as seen in Figure 16b. The 
time histories of pressures, Figure 18, show nonlinear behavior as the shock oscillates forward and aft with the pitch 
angle. 

Figure 18 is similar in format to Figure 17 and the reader is referred to the detailed description of Figure 17 to 
better understand the layout. The differences for the @ = 3.15° case are highlighted here. Snapshots of both the upper 
and lower surface pressure distributions are shown for the identified time points, color-coordinated with the symbols 
in the top graphic as discussed for Figure 17. Here, however, the upper surface snapshots are shown for all time points 
by the citrus colored region. 

For the cases at Q% = 3.15°, there is an upper surface shock, but the flow is not separated. Thus, the aftmost shock 
position occurs when the angle of attack is at the peak value. The movement of the shock is monotonic as the alpha 
value increases to that peak and monotonically decreasing as the angle of attack decreases away from that peak. This 
is shown by the snapshots as well as the color contour plot in the bottom graphic for each case. The three subplots 
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Figure 16. Time histories of analyzed conditions nearest the flutter conditions; Medium grid. 


show increasing dynamic pressure, with Figure 18 showing the behavior for an unstable case. The time history shows 
a growing amplitude oscillation and the snapshots show that the upper surface shock moves further forward with each 
oscillation through the negative aeroelastic angle range. Here, the total angle of attack is approaching 0°. On the 
contour plot, the supersonic region ahead of the shock is shown by the dark blue regions which transition abruptly 
to lighter blue regions as the pressure changes rapidly across the shock. The sinusoidal pattern formed by the edges 
of the dark and light blue regions shows that the pressures associated with the upper surface shock structure oscillate 
sinusoidally in time. 

The most thorough exploration of the dynamic solution results was done at a rigid angle of attack of 5°. Here, 
the flutter onset boundary was computed for fixed excitation amplitudes of 0.5 and 5.0. From the results shown in 
Figure 14, at gvelyp = 5, the system was unstable at g = 35 psf for the medium grid results and Figure 15c indicates 
that the initial positive direction oscillation angle of attack reached 6.7°. The phase relationship between the angle of 
attack and the shock movement shifted by 180° (i.e., there is a sign change). Linear interpolation between results at ¢ 
= 25 and g = 35 psf gave an estimate of flutter onset at g = 30 psf. 

The behavior of the pressure distribution for different parts of the oscillation cycle are presented for & = 5° cases 
at three dynamic pressures in Figure 19. For the Q& = 5° cases, flow separation intervenes and the shock movement 
is not monotonic and the aftmost position of the shock doesn’t occur when the angle of attack is largest. The aftmost 
shock location occurs when the total angle of attack reaches approximately 6.5°. This is observed most readily by 
examining Figure 19c, which shows the post-flutter behavior. The aftmost shock position is shown in the snapshot 
pressure distribution plot by the blue line. This snapshot corresponds to a time point near 201.35 seconds, identified 
by the blue square symbol in the graphics of the aeroelastic twist angle time histories. This time point is also indicated 
on the contour plot by the black dashed line. The black dashed line passes through the aftmost point on the dark blue 
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Figure 17. Pressure behavior during oscillations, a = 0°. 
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Figure 18. Pressure behavior during oscillations, o = 3.15°. 


region on the contour plot. As described previously, the dark blue regions correspond to the supersonic plateau ahead 
of the shock. As the angle of attack continues to increase during the oscillation, the shock movement changes direction 
and moves forward. The forward movement halts as the pitching oscillation cycle progresses; the direction of shock 
travel changes when the total angle of attack decreases back to below 6.5°. To better illustrate the movement of the 
shock with the angle of attack change, a three-dimensional plot for an example case is shown in Figure 20 for & =5°, 


qG = 100 psf. 


Frequency analysis of time histories of the pressures and the pitching moment coefficient show dynamics at the 
flutter frequency (or the aeroelastic primary coupled frequency) but also at two, three and five times the flutter fre- 
quency. Rather than being superharmonics due to data processing, features of the time history data can be traced to 


each of these frequency multiples as the shock reverses over portions of the oscillation cycle. 
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The bottom of the transonic dip of the flutter boundary appears to coincides with separating/reattaching flow 
physics. At Mach 0.74, using the medium grid, this occurs at Qa; = 6.5°, where the flutter simulations without flux 
limiter and with standard perturbation velocity induce flutter for every dynamic pressure, from 0 psf to 300 psf. This 
is not to say that the predicted flutter onset at 0 psf is physically meaningful. The nonlinear responses in the pitching 
moment coefficient indicate that the simulated flow fields are different from the flow fields predicted at the lower 
angles of attack. 


I. Matched Point & Single Degree of Freedom Analyses 


The flutter summary plots for the medium grid simulations, shown in Figure 14, include two additional cases that 
haven’t been discussed: matched point analysis and pitch-only analysis simulations at % = 5°, g =50 psf with an 
initial generalized velocity perturbation of 5.0. Simulating a matched point condition means having Mach number, 
dynamic pressure and Reynolds number all corresponding to the same physical condition. This analysis result, shown 
by the lavender-filled square on each subplot in Figure 14, indicates that there is little influence on the damping or 
frequency value as it closely matches the result for the nonmatched condition. 

The second additional case is shown in Figure 14 by the dark red squares. This simulation was performed using 
only the pitch degree of freedom in the structural model. The results, again, show little influence on the damping or 
frequency. This result shows that, for this specific case, the flutter mechanism is a single degree of freedom—pitch- 
only—instability. It has not yet been evaluated whether this is true for all of the cases, or if possibly this is only true 
for those cases that are seen to deviate from the upper arch of results shown on the damping plot, Figure 14b. If the 
latter is true, this would indicate a change in instability mechanism due to amplitude of perturbation, as well as angle 
of attack. 
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Figure 19. Pressure behavior during oscillations, % = 5°. 


2. Spatial Resolution Effects 


Dynamic aeroelastic analyses were also performed using the extra fine grid. Summary flutter results are shown in 
Figure 21 for three rigid angles of attack (5°, 6.5° and 8°) and two values of initial generalized velocity (0.5 and 5.0). 
The system was unstable at g = 75 psf; linear interpolation using results at g = 50 and g = 75 psf gave an estimate of 
flutter onset at g = 60 psf. This is substantially higher than the flutter prediction from the medium grid results, g = 
30 psf. The discrepancy between the two grid results is in contrast to the results produced for AePW-2 at Mach 0.74, 
O% = 0°, where the grid size had almost no impact on the predicted flutter onset condition. The result does agree well, 
however, with the trends observed in the AePW-2 results at Mach 0.85, & =5°, where the spatial resolution (grid size) 
has been shown to have a substantial effect on the flutter prediction.” The requirement for increased spatial resolution 
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Figure 20. Pressure behavior during oscillations, 0% = 5°, g = 100 psf. 


appears to be driven by the presence of separated flow. Additionally, at @ = 50 psf, where the system analyzed with 
the extra fine grid was stable, the positive direction oscillation angle reached a maximum of 7.02°. This indicates that 
the critical angle of attack is slightly higher than the critical angle predicted using the medium grid (@ = 6.5°). No 
experimental data set is available to verify either prediction, however. 
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Figure 21. Flutter results; Extra fine grid. 


3. Initial Perturbation Effect 


The typical flutter calculation in FUN3D begins with an initial excitation “kick” at the start of the dynamic aeroelastic 
simulation. In the current paper, this kick was applied in the form of an initial value of generalized velocity. The 
magnitude of the kick is somewhat arbitrary; however, it cannot be so large that it causes large grid deformation and 
possible solution failure. A large kick in the pitch degree of freedom may also change the aerodynamic state of the 
flow, leading to incorrect flutter prediction for systems assumed to be subject to small perturbations in freestream 
conditions. On the other hand, a kick that is too small will result in the solution taking a very long time to develop. 
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In reference 2, the notion that an initial excitation kick might not be required to obtain a flutter solution was explored. 
Rather than applying a fixed magnitude perturbation, reference 2 explored the concept that the change in going from 
the rigid steady solution to the coupled aeroelastic solution would be sufficient to perturb the system and assess flutter. 
The damping value for the case examined using this modified perturbation method was similar to the value extracted 
using the nominal perturbation (gvelp = 5) applied in the current work, implying that the step input associated with 
applying the aeroelastic coupling is approximately equivalent to a perturbation level of 5. 

In the current work, the influence of the initial perturbation level on the flutter onset boundary was assessed at 
O% = 0° and & =S° angles of attack at fixed dynamic pressures, using the medium grid. The damping values extracted 
for these cases are presented in Figure 23. At & = 0°, dynamic pressures of 152 and 169 psf were examined. The 
damping values for these cases are shown to be gradual, monotonic and nearly constant slope functions of the initial 
perturbation level. 

The flutter onset boundaries assessed in terms of the perturbation size at 5° angle of attack show more complicated 
dependence. As shown in Figure 22, time histories resulting from an initial assessment at a dynamic pressure of 100 
psf show significantly different characteristics for cases with values of gvelo = 0.5 and 5. The figure shows the smooth 
sinusoidal variation of the pitch angle, overplotted with the pitching moment coefficient. For the large perturbation 
level, the pitching moment coefficient shows nonlinear characteristics, particularly at the peak of the pitching moment 
coefficient, corresponding to the time points just after the angle has reached the maximum value. 
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Figure 22. Time histories 0% = 5°, 100 psf; Medium grid. 


To further explore the effect, additional levels of perturbation were evaluated. Values of gvelo = 0.5, 1, 1.5, 2, 3 
and 5 were used to excite the simulation; the damping results are presented in Figure 23. For gvelo greater than or 
equal to 1.5, the initial positive direction oscillation went to & = 6.5°, and the system destabilized. The lower two 
perturbation level cases (0.5 and 1.0) were stable. The damping values indicate that there is a distinct transition in 
the physics between the simulations at 1.0 and 1.5. As shown in Figure 10 for the @& = 5°, g=100 psf case, the static 
aeroelastic twist angle is 1.05°. A plot showing the absolute value of the dynamic portion of the aeroelastic twist angle 
is given in Figure 24. The orange line shows the results for the stable gvelg = 1 case, which has a first oscillation peak 
amplitude of 0.3°; the brown line shows the unstable results for the gvelp = 1.5 case, which has a first oscillation peak 
amplitude of 0.5°. Combining the rigid angle of attack (5°), the static aeroelastic twist angle (1.05°) and the dynamic 
peak angles for the two cases on either side of the stability boundary produces values of Qjorai = 6.35° and 6.55°. Note 
that the first is stable and the second is unstable, reinforcing the interpretation that cases with a maximum angle of 
attack less than 6.5° are stable, while cases that exceed this angle are unstable. 


4. &=6.5° Cases 


For Q& = 6.5°, out of all the simulations performed, there is only one solution that might be considered as neutrally 
stable. All the others are unstable. Two example results from the medium grid simulations, both showing erratic 
progression and instability, are shown in Figure 25. The extra fine grid at g = 25 psf with a gvely of 0.5 shows 
an oscillation in pitch angle that, while oscillating erratically, doesn’t diverge over the cycles performed. Results 
discussed earlier for the extra fine grid also indicate that the critical angle is in excess of 7°, so this case may not have 
sufficient excitation to achieve the critical angle. In general all of the simulations at @& = 6.5° contained long and 
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Figure 23. Flutter results; Medium grid; Damping variation with initial perturbation level. 
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Figure 24. Flutter Results, a = 5°, 100 psf; Medium Grid; Variation of the absolute value of the dynamic portion of the 
aeroelastic twist angle with initial perturbation level. 


variable transient regions; many of the solutions also appear to change macroscopically due to computer round-off. 
These are matters for continued exploration. 


V. Continuing Analyses 


The companion publication to this work, reference 2, indicates that the back side of the flutter bucket (transonic 
dip) for the 5° angle of attack condition is near Mach 0.82. The flutter boundary across the transonic Mach range at 5° 
angle of attack, produced in reference 2, is reproduced here in Figure 26. As the angle of attack increases, the transonic 
dip of the flutter boundary has historically been observed to shift to the left (i.e., to a lower Mach number), exemplified 
by Figure 1. Based on the steady rigid results, 8° was investigated as a candidate angle of attack to see if the transonic 
dip would be shifted to below Mach 0.74. That is, & = 8° was explored to see if the flow was sufficiently separated 
such that the transonic dip of the flutter boundary would be moved to a lower Mach number. If so, this would produce 
a higher dynamic pressure for flutter onset, as was observed at Mach 0.82 for Q& = 5°. The extra fine grid was used 
for this investigation. The results did not show the desired shift to a lower Mach number, with the damping values 
indicating that the 8° system is not more stable than the 6.5° system for corresponding dynamic pressures. Thus, the 
back side of the flutter bucket is predicted to be at a Mach number greater than 0.74 for Q& = 8°. Higher angles of 
attack could be examined, but it seems likely that we have already applied the URANS analysis to conditions beyond 
the range where they are credible. Without experimental data to benchmark against, further exploration extending the 
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Figure 25. Example time histories, 0 = 6.5°; Medium grid. 


angle of attack range above 8° doesn’t seem fruitful. 
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Figure 26. Flutter boundary assessment, Mach 0.6 - 0.85, a = 5° using FUN3D, from reference 2. 


The results shown in Figure 26 overlap the results in the current paper at Mach 0.74. Note that in these results, 
the back side of the transonic dip for & = 5° is predicted to begin by Mach 0.82, as indicated by the increased 
flutter dynamic pressure. The small triangular symbols in Figure 26 show the experimental points that were produced 
at nearby angles of attack. The measured angles of attack for these unstable conditions were 5.3°, 5.4° and 5.5°. 
No mention was made in the publications of the experimental data that a flutter-free pass was made below these 
conditions. This suggests that the flutter onset condition could be even lower than the instability points reported. 
These experimental points give some confidence that the flutter predictions in this range are reasonable, although this 
comparison is far from conclusive. 

Flutter computations are currently being performed at higher angles of attack at Mach 0.74 to attempt to identify 
the edge of the transonic flutter bucket, that is, the point at which the flow remains separated and the flutter onset 
condition rises. At Mach 0.80, simulations are being performed at angles of attack between 5° and 10° to attempt to 
establish the onset of the back side of the flutter bucket. Additional studies that may be pursued include modeling 
the system with a Lattice Boltzman approach, explicit modeling of the splitter plate, more detailed examination of the 
influence of the tip vortex, matched point condition effects, mapping the flutter envelope explicitly using only the pitch 
degree of freedom, and revisiting the two-dimensional airfoil configuration to better understand the influences of the 
wing tip vortex. 
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VI. Conclusions 


The transonic dip of the flutter boundary encountered by the BSCW is due to separating/reattaching flow. Using 
a medium grid with 9 million nodes, the critical angle for this condition at Mach 0.74 is predicted to be 6.5°. The 
predicted angle is slightly larger using an extra fine grid with 35 million nodes. Effects of this critical angle are seen 
in the rigid steady results as a reversal in the shock position with increasing angle of attack. The effects are seen in 
the static aeroelastic solutions as reductions in the total angle of attack with increasing dynamic pressure when the 
total angle of attack exceeds 6.5°. The effects seen in the dynamic aeroelastic solutions are a change in the phase of 
the FRFs of the pitching moment coefficient when the dynamic pressure exceeds the flutter condition. Examination 
of the shock movement over a cycle showed the same reversal as observed in the steady rigid and static aeroelastic 
solutions. This behavior leads to different frequencies being significantly represented in the responses. This isn’t a 
superharmonic effect, but rather the effect of separation onset part way through the cycle. 

Different perturbation levels used to provide initial excitation at the start of the dynamic analyses give different 
flutter boundaries when separation onset is part of the flutter inducement (i.e., within the transonic dip region). When 
separation is not a factor, the perturbation level matters little. At Mach 0.74, the transonic dip of the flutter boundary 
is defined by solutions perturbed such that they reach 6.5° total angle of attack. At other transonic Mach numbers, it 
is anticipated that there are similar critical angles. It is also anticipated that the critical angle of attack will be lower at 
higher Mach numbers and higher at lower Mach numbers. Additional experimental quantitative information is needed 
to assess where the simulation represents physical reality. 
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